Contact process on a Voronoi triangulation 
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Abstract 

We study the continuous absorbing-state phase transition in the contact process on the Voronoi- 
Delaunay lattice. The Voronoi construction is a natural way to introduce quenched coordination 
disorder in lattice models. We simulate the disordered system using the quasistationary simulation 
method and determine its critical exponents and moment ratios. Our results suggest that the 
critical behavior of the disordered system is unchanged with respect to that on a regular lattice, 
i.e., that of directed percolation. 
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I. INTRODUCTION 



Nonequilibrium phase transitions between an active (fluctuating) state and an inactive, 
absorbing state arise frequently in interacting particle m ode.s g chemical cata.ysis Q, 
interface growth [3J , epidemics [4J] and related fields. In spatially extended systems, exem- 
plified by the contact process such transitions are currently of great interest, which has 
been heightened by recent experimental confirmations of absorbing-state phase transitions 
in a liquid crystal system [6], and in a sheared colloidal suspension |7']. Much of this work 
is focused on issues of universality, aimed at characterizing the critical behavior of these 
models in terms of universality classes , Q. It has been conjectured 
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12j that 



models with a positive one-component order parameter, short-range interactions, and ab- 
sence of additional symmetries or quenched disorder belong generically to the universality 
class of directed percolation (DP), which is considered the most robust universality class of 
transitions to an absorbing state. 

The contact process (CP) is one of the simplest and most studied models belonging to 
the DP universality class. Of particular interest is how spatially quenched disorder affects 
its critical behavior [lsj]. Quenched disorder, in the form of impurities and defects, plays 
an important role in real systems, and may be responsible for the rarity of experimental 
realizations of DP [14]. Quenched disorder in the contact process on a regular lattice has 



been studied in the forms of random deletion of sites or bonds 
spatial variation of the control parameter 
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16l. Il7i , and of random 



20] . All these studies report a change in 



the critical behavior of the model. These findings are consistent with Harris' criterion [21], 
which states that quenched disorder is a relevant perturbation if 

du ± < 2, (1) 

where d is the dimensionality and u± is the correlation length exponent of the pure model (In 
DP this inequality is satisfied in all dimensions d < 4, since u± = 1.096854(4), 0.734(4) and 
0.581(5), for d — 1, 2 and 3, respectively 22j, |23|, |24|.) Some controversy remains whether 
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251 ]. or whether they change 



the exponents change continuously with degree of disorder 
abruptly to the values in the strong disorder limit corresponding to the universality class of 
the random transverse Ising model, as suggested by Vojta in a recent work 19]. 



Harris' criterion determines the relevance of disorder in the form of independent random 
dilution (of sites and/or bonds) in a regular lattice. A somewhat different situation arises 



when the underlying graph is not periodic, as is the case in a deterministic aperiodic struc- 
ture, or in a graph with a random neighbor structure such as the Voronoi triangulation. 
To determine the relevance of disorder in these cases, the following heuristic extension of 



Harris' criterion was proposed by Luck 26] : Consider a spherical patch Q with radius R on 
a given realization of a graph. The patch encloses a number B(R) of vertices, which scales 
as B(R) ~ R d . The average coordination number in the patch is given by 

Let the fluctuation of the coordination number around its expected value, J = q, decays as 
^(j) = (\ J (R) -Jo\) ^ {B{R)) -[i-») „ R -*o.-») i (3) 

Jo 

when R — > oo. Here, u is defined as the wandering exponent of the triangulation. Nearby 
the critical point A = (A — A c )/A c = 0, the fluctuations cr^(J) of the average coordination 
number in a correlation volume scale as 

ot(J) ~ CI d/2 ~ A" ±d/2 , (4) 

since £j_ ~ A~ I/± . Considering a large correlation volume, R ~ the resulting shift of 
the critical point, induced by the fluctuations in a correlation volume is proportional to 



^du±(i-(j) y/ vai (q.^ where var(gj) = (qf) — (qi) 2 . Then, in order that the regular critical 
behavior remain unchanged, these fluctuations should die out when A — >• 0, which is true if 
us does not exceed a threshold value given by 

Thus, in principle, the Harris-Luck criterion permits one to predict the effects of quenched 
disorder in models defined on structures such as quasi-crystals or even random lattices. 
(Note that for independent dilution, u> = 1/2, and Luck's expression reduces to the Harris 
criterion.) 

In this work we investigate whether disorder in the form of a quenched Poissonian coor- 
dination disorder alters the critical behavior of the contact process, by studying the critical 



behavior of the process on a Voronoi-Delaunay (VD) type random lattice 271, |28[]. The VD 



lattice represents a natural way of introducing quenched coordination disorder in a lattice 



model, and also plays an important role in the description of idealized statistical geome- 
tries such as planar cellular structures, soap throats, etc. 27|. In this lattice, the sites are 



spatially distributed following a Poisson distribution, and the coordination number q varies 
randomly, with 3 < q < oo and q = 6 in the infinite-size limit. Our results suggest that 
coordination disorder does not change the critical behavior of the contact process. 

The balance of this paper is organized as follows. In the next section we review the 
definition of the contact process and detail construction of the VD lattices as well the 
simulation methods used. In Sec. Ill we present our results and discussion; Sec. IV is 
devoted to our conclusions. 



II. MODEL AND METHOD 

Consider a bounded domain Q in a d- dimensional space in which N nodes are randomly 
placed with uniform distribution. The Voronoi diagram of this set is a sub-division of the 
domain into regions Vi (with i = 1,2,..., N), such that any point in Vi is closer to node i 
than to any other node in the set. Figure 1 (a) shows a patch of a Voronoi diagram. The 
points whose cells share an edge are considered neighbors. The dual lattice, obtained by 
linking neighboring sites is the Voronoi-Delaunay network, exemplified in Fig.l (b). One of 
the characteristics of the dual lattice is that its local coordination number varies randomly, 
with the distribution shown in Fig.2. In this work, we take periodic boundary conditions, i.e, 
the domain Q has a toroidal topology. In order to construct the lattices we follow the method 
of Ref . 29( . For simplicity, we express the length L of the domain Q in terms of the size of 
a regular lattice L = \/~N. 

The CP, deduced a S a , oy m odeP for epide^c S p r ead lng fl , is a stochastic 

interacting particle system defined on a lattice, with each site either occupied (&i(t) = 1), or 
vacant (<Ji(t) = 0). Transitions from <7j = 1 to <7j = occur at a rate of unity, independent of 
the neighboring sites. The reverse transition is only possible if at least one of its neighbors 
is occupied: the transition from <jj = to cr, = 1 occurs at rate Ar, where r is the fraction 
of nearest neighbors of site i that are occupied; thus the state a, = for all i is absorbing. 
(A is a control parameter governing the rate of spread of activity.) 

In the simulation we employ the usual simulation scheme jl| , in which annihilation events 
are chosen with probability 1/(1 + A) and creation with probability A/(l + A). In order to 



FIG. 1: (a) A patch of a Voronoi Diagram, (b) The corresponding dual lattice to the diagram 
shown in (a). (Color online). 
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FIG. 2: Degree distribution P(q) of the Voronoi-Delaunay lattice, for system size L = 2560. 

improve efficiency, the sites are chosen from a list of currently occupied sites. In the case 
of annihilation, the chosen site is vacated, while, for creation events, one of its q nearest- 
neighbor sites is selected at random and, if it is currently vacant, it becomes occupied. The 
time increment associated with each such event is At = l/N occ , where N occ is the number of 
occupied sites just prior to the attempted transition. 

In the studies reported here we sample the quasistationary (QS) distribution of the pro- 
cess, (that is, conditioned on survival), which has proven a very usefu. too. in the study 



of processes with an absorbing state hj, l30|, l3l| . For this purpose, we employ a simulation 



32] 



method that yields quasistationary (QS) properties directly, the QS simulation method 
The method is based in maintaining, and gradually updating, a set of configurations visited 
during the evolution; when a transition to the absorbing state is imminent the system is 



instead placed in one of the saved configurations. Otherwise the evolution is exactly that of 
a conventional simulation. 



III. RESULTS AND DISCUSSION 



We performed extensive simulations of the CP on Voronoi-Delaunay random lattices of 
L = 20, 40, 640, using the QS simulation method. Each realization of the process is 
initialized with all sites occupied, and runs for at least 10 8 time steps. Averages are taken in 
the QS regime, after discarding an initial transient which depends on the system size. This 
procedure is repeated for each realization of disorder (For each size studied, we performed 
averages over 200-300 different lattices). 
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FIG. 3: Quasistationary density of active sites p as a function of the control parameter A. System 
sizes: L = 20, 40, 80 and 160, from top to bottom. 



In Fig. 3 we show the quasistationary density p as a function of the control parameter A 
for several values of L. We observe, as expected, a continuous phase transition from an active 
to an absorbing state. Since, due to topological constraints, the Voronoi-Delaunay lattice 
posses q ~ 6, the value of the critical point is shifted from A c = 1.64877(3) 33| (regular 
square lattice) to A c = 1.54266(4) (the increase in q facilitates creation). This is very close to 
the critical value, A c = 1.54780(5), for the regular triangle lattice, obtained using the same 
methods as described below. It is notable that the critical value of the disordered system 
is about 0.3% smaller than that of the regular lattice with the same average connectivity. 
Fig. 4 shows how the QS density of active sites varies with the coordination number q. 
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FIG. 4: QS density of active sites p(q) versus A, for sites with q = 3, 4, ...,10, from bottom to 
top. Circles: average over all sites. System size L = 160. 

At the critical point we find that the quasistationary density decays as a power law, 
p ~ L~^ x , as shown in Fig. 5. Our simulation data follow a power law with the exponent 
fl/v± = 0.791(7), while the value for DP in two spatial dimensions is 0.797(3) 331 ] . 
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FIG. 5: QS order parameter p versus system size L at criticality (A = 1.54266). 



Another important quantity is the lifetime of the QS state, r. In QS simulations we 
take r to be the mean time between successive attempts to visit to the absorbing state. 
Fig. 6 shows that at the critical point, the lifetime also follows a power-law, r oc L z , with 
z = vu/ux = 1.78(3), as compared with the literature value of 1.7674(6) for the DP class 



33]. 
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FIG. 6: Critical lifetime of the QS state r versus L. 

Complete characterization of a nonequilibrium universality class requires the determina- 
tion of at least three independent critical exponents. To this end we perform initial decay 
studies on large systems, starting with a fully occupied lattice. While the CP with random 
dilution exhibits a logarithmic relaxation [16|], on the VD lattice we observe a clear power- 
law decay. Finite size scaling in this case predicts that p ~ t~ s . A least-squares fit for the 
data shown in Fig. 7 yields S = 0.453(9), in very good agreement with the standard value 
of 5 = 0.4523(10) for DP |33|. 




FIG. 7: Decay of the order parameter starting from a full lattice of size L = 2560. A = 1.5428, 
1.5427 and 1.5426, from top to bottom. (Color online). 



Moment ratios (or reduced cumulants) represent an alternative method for identifying 
the universality class of a continuous phase transition 3J, |35|, |36j . Here we analyze the 



critical moment ratio m = (p 2 )/(p) 2 - This quantity is analogous to Binder's reduced fourth 



cumulant 



371 ]. at an equilibrium critical point: the curves m(X,L) for various L cross near 



A c (the crossings approach A c ), so that m assumes a universal value m c at the critical point, 
as can be seen in Fig. 8. In this case, our data yield a universal value of m c = 1.328(6), again 
in very good agreement with the best known value for the CP on a regular square lattice, 
m c = 1.3257(5) |34j. 




FIG. 8: Quasistationary moment ratio m versus InL, for A = 1.54256, A = 1.54260, A = 1.54264, A = 
1.54268 and A = 1.54280, from top to bottom. System size: L = 640. Inset: Quasistationary 
moment ratio m versus p, system sizes: L = 20,40,80, 160. (Color online). 



In summary, our results reveal that the absorbing phase transition of the contact process 
defined on Voronoi-Delaunay random lattice belongs to the directed percolation universality 
class. These results are somewhat surprising, since the wandering exponent for these lattices 
was numerically evaluated in a extensive work by Janke and Weigel 38j, who found that 
lo = 1/2, i.e, the relevance criterion for such lattices reduces to the usual Harris criterion, 
eq.(l). 

In the equilibrium context the Harris-Luck criterion has been verified numerically on 



random latices in several models, such as the Ising model |39l . |40| and percolation 41 ] 



However, Monte Carlo simulations for the q = 3 Potts model [42j as well for the Ising model 



in 3D 



43 



44j and for the spin-3/2 Blume-Capel model 45] yield results that contradict the 



relevance threshold given by the Harris-Luck criterion. Simulation results for nonequilibrium 
models, viz. the majority- vote model on a random lattice 46j, also appear to contradict 
this relevance criterion. 
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FIG. 9: Survival probability versus time, in the critical CP on a Voronoi lattice (solid line), and for 
the critical CP on a square lattice with random dilution of 2% (doted) and 5% (dashed). System 
sizes: L = 640 (left curves) and L = 1280 (right). 



In Refs. [38|, |44j it is suggested that Voronoi disorder appears not to alter the critical 
behavior because it is intrinsically weak, and that the usual hallmarks of quenched disorder 
would in fact manifest themselves in larger systems. In order to test this hypothesis, we 
compare in Fig. 9 the survival probability Ps (starting with a fully occupied lattice) of the 



CP on the VD lattice and on a regular (square) 



attice with weak dilution. It is known 



that the diluted CP exhibits activated disorder [19J, due to emergence of favorable regions, 
leading to logarithmically slow dynamics [if]]. We find that while in the CP on the VD 
lattice the survival probability decays exponentially (as in the ordinary contact process), in 
the diluted CP the behavior is clearly different. 

Notice that the effect of the "rare regions" is clearly visible for the system sizes used 
here, even for the smallest dilution (2%): the decay of the survival probability is clearly 
slower than exponential. On the square lattice with dilution i < 1, the variance of the 
connectivity var(q) ~ Ax, so that var(q) ~ 0.08 for dilution 0.02. This is less than 5% 
of that for the VD lattice, where var(q) = 1.779(2). We also should mention that the 
effects of quenched disorder in CP are stronger than in the three-dimensional Ising model: 
in the latter, quenched disorder provokes a difference in the second digit in the exponent v 



471 ]. while for models in the DP class even weak disorder changes the critical dynamically 



properties drastically 



48|] 



IV. CONCLUSIONS 

We performed large-scale simulations of the contact process on a Voronoi-Delaunay ran- 
dom lattice, which exhibits quenched connectivity disorder in the model. Our results suggest 
that this kind of disorder does not alter the DP character of the transition, in contradiction 
with the Harris-Luck criterion. Given the large systems sizes and long simulation times 
used, it appears unlikely that the system will eventually cross over to non-DP scaling. Thus 
it remains an open question why an argument of the Harris-Luck type is not applicable in 
some cases. Our results also reveal that the DP universality class may be even more robust 
than asserted in the usual DP conjecture, in the sense that not all kinds of quenched disorder 
are relevant perturbations. 
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